Adaptive Sampling Approach to the Negative Sign Problem 
in the Auxiliary Field Quantum Monte Carlo Method 
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We propose a new sampling method to calculate the ground state of interacting quantum systems. 
This method, which we call the adaptive sampling quantum monte carlo (ASQMC) method utilises 
information from the high temperature density matrix derived from the monte carlo steps. With 
the ASQMC method, the negative sign ratio is greatly reduced and it becomes zero in the limit At 
goes to zero even without imposing any constraint such like the constraint path (CP) condition. 
Comparisons with numerical results obtained by using other methods are made and we find the 
ASQMC method gives accurate results over wide regions of physical parameters values. 
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The negative sign problem in quantum monte carlo 
simulations has been an extremely serious problem. In 
the case of the auxiliary field quantum monte carlo 
(AFQMC) method applied to the two-dimensional 

(2D) Hubbard model, we always have this problem in 
the ground state and in low temperature regions except 
at half-filling, where we have particle-hole symmetry. 
The negative sign difficulty is somewhat reduced and cal- 
culations become more feasible only when the filling is 
far away from half- filling, and/or when we use a smaller 
value of U/t, and/or when we use a smaller value of the 
inverse temperature P or the projecting time r. Such 
a reduction is also observed when the electron filling is 
such that the corresponding non-interacting model has 
a closed-shell electronic structure and the system size is 
small. The quantum monte carlo results on the 2D mod- 
els and the Hubbard ladder model obtained so 
far are in regimes where the above conditions are satis- 
fied. However, to get more insight into the physics of 
strongly correlated electrons, it is highly desirable to de- 
velop more robust numerical methods whose applicability 
and accuracy do not depend on the details of the system 
to be studied. 

Recently, Zhang, Carlson and Gubernatis developed 
the constrained path quantum monte carlo (CPQMC) 
method. |l^ Their method consists of the two ideas: the 
random walker in the configuration space (RWCS) and 
the constrained path (CP) conditions, which is local. The 
latter is a variant of the nonlocal positive projection con- 
dition proposed by Fahy and Hamman but is simpler 
to impose. Both of their methods are variational but the 
rate of convergence to the ground state of the CPQMC 
method against a variation of the trial wavefunction has 
been reported to be very fast in weak and intermediate 
regions of U/t and the closed shell case. We will demon- 
strate that the CP condition is not necessary to reduce 
the negative sign ratio, if we adopt the adaptive sampling 
method which we propose here in the standard projector 



AFQMC (PAFQMC) algorithm. 

In applying the PAFQMC method to the Hubbard 
modd: H = T + V2,T ^ -tE(.,,>(4c,. + i^-C), 
V2 = U J2i ^iiT^iii where {i,j) indicates the sum is taken 
over the pairs of nearest-neighbor sites, we use the den- 
sity matrix: (^ftl exp (— tTJ) IvPt). The density matrix is 
expressed by using the Suzuki- Trotter (ST) formula and 
the Stratonovitch-Hubbard (SH) transformation in the 
following way: 

(*t|cxp {~TH)\^t) = Tr,(i),,(2),..,,(L)(*nl^T('^(i)) • • • 
ST(a(l))|**T) X (*u|Si(a(L))...Bi(a(l))|vI/,^), 

where L — r/Ar and At is the discretized project- 
ing time r, Ba{l) — exp (— ArT) exp (— AtViq(O) 
Vic,{l) = 5,j{aaucr,{l) - (l/2)Art/). In the Vi^il) factor 
ajj is defined to be tanh^^ tanh{ATU /A) and a = +1 
for up spins and —1 for down spins. Hereafter, we de- 
note [/^(r, 0) — B{cr{L)) ■ ■ ■ B{a{l)) and the spin in- 
dices will be suppressed. Products of both up-spin and 
down-spin elements are implicit throughout this letter. 
T is the total projection time. The trace over the Ising 
SH fields: a{l),a{2) ■ ■ ■ ,a{L) is achieved by using the 
importance sampling method. The site index i of the 
SH fields is suppressed here. After we take the trace 
over the SH fields, we can calculate any ground state 
expectation value and the ground state wavefunction: 



t) 



Tr. 



tT(l),<T(2),...,<T(L) 



B{a{L))---B{a{imn)- 



There remains freedom how we take the trace. In the 
standard AFQMC method, the trace over the SH fields 
is taken simultaneously and we use the weight func- 
tion: P — |(5't|[/o-(T, 0)|\['t)|. There remains another 
choice in how we take the trace. If we define: = 
rr,(,)B(a(z))|^',_i), and 4-0 = 

\^ exact) ^ Tr,(2),-ML)B{<T{L)) ■ ■ ■ B(a(2))|*i) 
= rr,(3),...,.(L)SKi))---i?(a(3))|vI/2) 
= ---=Tr„^L)B{c7{L))\^L-i). 
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We can then take the trace sequentially. This is similar to 
the propagation process in RWCS and both this method 
and the standard AFQMC method should give identical 
results, if proper monte carlo sampling is done. 

To take the sequential trace properly, we introduce the 
adaptive sampling method and we define the weight func- 
tion as follows: 

Pr' = |(*t|C/<.(T,T')|*t)(*t|C/a(r',0)|*t)|,(0 <t'< t), 

rather than the standard weight function: P = 

K^tlUrrir, 0)|4't)|- t' is an imaginary time whose SH field 
is tried to be updated and r' changes during the simula- 
tion; the SH field to be updated belongs to the I = t'/At- 
th imaginary-time slice. We use the local-flip update 
scheme and the field is updated from I = 1 to I = L 
( from t' = At to r' = r). The ratio used in judging 
to accept or to reject the update of the SH field of the 
i-th site and of the I = t' /At -th imaginary-time slice is 
given as follows: 

r.. = |(vI/,|C/_,^(,,)(r',0)|vI/,)/(VI/,|t/+^^(,,)(r',0)|^'i)|. 

We use the heat bath method and the acceptance prob- 
ability is given by Rr' = rr' /{I + rr')- The r' dependent 
weight function is the unique feature of our adaptive sam- 
pling method. The update of the SH field is initiated by 
using the "high temperature density matrix" (the pro- 
jecting time r' ~ At ). All the SH fields on the 
imaginary-time r' are tried to be updated, where A^ is 
the number of sites. The series of trial updates on the 
imaginary-time r' is repeated M times. After the tri- 
als have completed, r' increases by At. The procedure 
is repeated until r' becomes r. All these constitute a 
sweep. The next sweep starts with t' = At, again. Wc 
call our method the adaptive sampling quantum monte 
carlo (ASQMC) method. 

Measurements may be done with the re-weighting func- 
tiom Rw{t') = \WJW,{t% = {^t\U^{TM^t), 
W,{t') = (*t|[/,(T,T')|*t)(*t|?7.(r',0)|*t). Let O be 
one-body operators such like: O ~ ct c^^, where k is the 
wave vector and s is the spin variable. The expectation 
value of O can be calculated as follows: 



Tr,{0)SignW„Rw{T'){\W,{T')\/Tr,\W„{T')\) 
Tr^SignW^Rw{T')QW^{T')\/Tr^\W^iT')\) 

We may make measurements only by use of configura- 
tions of the SH field obtained just after the trial updates 
at t' = T (Z = L) and then the re-weighting function 
Rw{t) = 1. In this case, we may evaluate the expecta- 
tion value as follows: 



where 



^ ^ (*t|C/.(T,0)|vI/i) 



and Us is the number of samples. For the two-body (or 
more) physical quantities, we use the Wick theorem. The 
measurements in this form are called a mixed estima- 
tor: (*t|0|*e2;act)- Onc of thc difficulties of the ASQMC 
method is measurements because the mixed estimator 
does not give exact expectation values other than for en- 
ergy and sign. To remedy the situation, we may use the 
following weight function: 

P=\{^t\Ua{T,0)\^t)\{0<T' <Tc), 



Pr' = |(^t|C/.(T,T')|*t)(*t|t^<x(T',0)|*t)|, (Te < t' < t). 

Measurements in this case are made only when < t' < 
Tc . In this interval Rw{t') = 1 and measurements such 
as: 



where 



{0) = 



{^t\Ua{r,T')OU,{T', 0)\^t) 
(vI/ilf/^T.O)!*^) 



may be done there. We call this measurement the stan- 
dard measurement. 

Instead of using ('''') to define Rw{t'), we may be 
able to use any weight function which docs not necessary 
have any physical correspondence to the original prob- 
lem. For example, we may use the weight function of the 
half-filled Hubbard model fl^ to study the doped Hub- 
bard model. The re-weighting function in this case is 
Rw = \Wa/^cr\- The expectation value of any one-body 
operator is given: 

m) = 

Tr,{0)SignW^Rw{\n„\/Tr^\n„\) 
Tra SignWaRw /Tr„ \n^\) 

As the re- weighting procedures do not introduce any ap- 
proximation to the theory, both of the mathematical ex- 
pressions for the expectation value of operator O are 
exact. The numerical effectiveness of the re- weighting 
methods depends on whether SignW^ is mostly positive 
and if thc re-weighting function Rw is not very small. 
Unless the two conditions are satisfied, numerical appli- 
cation of re-weighting methods will not be successful. 

The ergodicity is guaranteed in our method. The 
ASQMC method is one of the exact algorithms like the 
standard AFQMC method and some other re-weighting 
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methods applied to it, irrespective of the sign prob- 
lem. This can be checked numerically at half-filling, 
where there is no sign problem and comparisons of re- 
sults obtained by using the ASQMC method and the 
standard AFQMC method are possible. We compared 
calculated results of the spin-spin correlation function: 
S{q) = l/NJ2^,Jel^p[-iq■ in -fj)]SfS^ of the 4x4 
Hubbard model at the half-filling. The results are Fourier 
transformed into the real space and they are plotted 
against the distance R in Fig. |l|. We find nice agree- 
ment over all the distances within the cluster. 

With our adaptive sampling method, we improve the 
trial wavefunction ^'^(i = 0, 1, • • ■ , L — 1) sequentially so 
that the final short time projection: exp(— ATiJ)|^i_i) 
extracts out the ground state and therefore ASQMC 
method is expected to be less affected by the negative 
sign problem. Our method is easily implemented by mod- 
ifying the standard AFQMC code. We do not use the 
population control procedure in our implementation so 
that At needs to be small. Here, we have studied the 
minimal case M = 1. 

The benefit of our ASQMC method is that it greatly 
reduces the negative sign ratio without imposing any con- 
straint such like the CP condition. We calculated the 
expectation value of the sign (Sign) as a function of At. 
Tc was set to be zero and we used the mixed estimator 
method to calculate the sign expectation value {Sign). 
We used the projecting time r/t = 10. Calculations were 
made with various physical parameters values of the 4 x 
4 Hubbard model, but we show the result obtained with 
U/t = 8 and with the density p = N^/Ns = 0.875, where 
there is very serious negative sign problem when we use 
the standard AFQMC method. The result is shown in 
Fig. 1^. When At — > 0, {Sign) — * 1. We have no negative 
sign problem in the limit of At — > 0. Over a wide range 
of finite At, the negative sign ratio obtained with our 
ASQMC method is much reduced in comparison with the 
standard AFQMC method. When we use smaller value of 
U jt or when we study the closed shell filling case, {Sign) 
tends to become much closer to 1 with the same value of 
At. The rate of convergence to {Sign) — > 1 depends on 
values of physical parameters. 

Our ASQMC method is expected to give accurate nu- 
merical results, because SignWa is mostly positive and 
Rw{t') = 1 . ( Measurements are done either at t' = t 
or < t' < Tc . ) We prefer to use the standard measure- 
ment here to reduce the At error, which comes not only 
from the Trotter error but also from the finite step size 
At used in taking the sequential trace of the auxiliary 
fields. The ground state wavefunction obtained by using 
finite At may be expanded by powers of At: ^'(At) ~ 
■^,^act + At*'(0), where *'(0) = (a/aAT)*(AT)|Ar=o 
and \l/'(0) is orthogonal to exact- If we use the 
mixed estimator measurement, we obtain the expression: 

(*t|F|*(AT)) = Eexact{-^t\^ exact) + At {-^ t\H \^ ' {Q)) . 

On the other hand if we use the standard measurement, 



we get : {^ {AT)\H\<i! {At)) = Eexact{^ exact I ^ exact I ~r 
At'^ {^' {{))\H\^' {O)) , where Eexact is the exact ground 
state energy of the Hamiltonian H. We notice that larger 
At error remains in the mixed estimator measurements 
than in the standard measurements, which is also ob- 
served in numerical simulations. So the standard mea- 
surement is the better choice even when one calculates 
energies. The negative sign ratio obtained by using the 
standard measurement method is somewhat larger than 
that obtained by using the mixed estimator method, but 
it is still far more reduced than that obtained by using 
the standard AFQMC algorithm. Hereafter throughout 
this article, we adopt the standard measurement method 
in our ASQMC calculations. 

We compared calculated total energies, one and 
two-body correlation functions obtained by using the 
ASQMC method with those obtained by using the Lanc- 
zos exact diagonalization (ED) method and other quan- 
tum monte carlo methods. Values of At used are 0.025 
and 0.0125. A 5000 sweep run on the 8x8 lattice takes 
almost 80 hours of CPU time on the Alpha workstation 
with the Alpha 21164 / 533MHz CPU chip. We found the 
ASQMC method gives accurate results over wide range 
of physical parameters values on 4 x 4, 6 x 6, and 8x8 
lattices. In the following paragraph, we show results of 
such comparisons. 

First, we compare total energies oi U/t — 4 Hubbard 
model of various lattice size: 4 x 4, 6 x 6, 8 x 8. Fillings 
p are such that both the closed shell case and open shell 
case are included. Agreements of the total energies calcu- 
lated by using the ASQMC method with those obtained 



by using ED, |T3| the CPQMC, |n|l and the AFQMC |] 
methods are rather nice up to 3 digits over the wide range 
of fillings and lattice size. This is shown in Fig. ^. To 
see this in more detail, we have plotted relative errors of 
the total energies in Fig. ^. The CPQMC results agrees 
with our ASQMC results within 0.1% of errors in the all 
cases studied here, but the AFQMC results does not al- 
ways agree with our ASQMC results very precisely. The 
largest error we found is 0.5%. It seems that the extrap- 
olation procedure used to make inferences of the ground 
state energies without taking in to account of statistical 
errors described in Ref. does not work so well. 

We next compare total energies of the 4x4 Hub- 
bard model for various U jt values calculated by using 
the ASQMC, the CPQMC ^ and ED methods. |l|] 
We set the filling p = 14/16. Agreements between the 
ASQMC and ED results are very good but we found 
small systematic increase of deviation of the CPQMC 
results from the other two results, as we increase U jt. 
To see in more detail this tendency found in Fig. |^, we 
have plotted relative errors of the total energies in Fig. ||. 
While the CPQMC results deviate from the ED results 
systematically, the ASQMC results do not. Because we 
do not impose any constraint such like the CP condition, 
the ASQMC method does not have any systematic bias 
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even in the large U /t region. This demonstrates that our 
ASQMC method is superior to the CPQMC method in 
the large U/t region. 

To provide other benchmark of our ASQMC method, 
we have calculated the d-wave pairing correlation func- 
tion of the 4x4 Hubbard model defined by: < 
0(i?)Ot(o) > with 0'^{R) = c\{R)c\(R + x) + 

c\{R)c\{R - x) - c\{R)c\{R + y)- c|(i?)c|(i? - y) and 
compared it with the result obtained with the Lanczos 
diagonalization method. The result is shown in Fig. 

. We again, obtained a good agreement with the Lanc- 
zos diagonalization result over all the distances within 
the cluster studied. 

As far as the present author knows, there are no other 
reliable numerical methods than ED, the AFQMC and 
the CPQMC methods to compare with our ASQMC 
methods for the 2D Hubbard model. Our ASQMC re- 
sults always agree well with the best of the results ob- 
tained by using the other methods. Thus our ASQMC 
method turns out to be a very accurate method. The 
ASQMC method will be useful to study a wider physical 
parameter region including the large U jt region, many 
more cases of filling and band structures, including the 
open shell cases, and larger lattice sizes that have not yet 
been explored by the quantum monte carlo methods. 

To summarize, we have proposed the adaptive sam- 
pling method to utilize information from the "high tem- 
perature density matrix" in the thermalization process 
of monte carlo steps to calculate the ground state. With 
the adaptive sampling method, the negative sign ratio 
decreases to when At — > without imposing any con- 
straint such as the CP condition. Over a wide range of 
finite At, the negative sign ratio is far more reduced than 
that in the standard AFQMC method. We compared 
calculated energies and two-body correlation functions 
obtained by using our method and with those obtained 
by using the Lanczos diagonalization method and other 
quantum monte carlo methods found in the literature and 
we found the ASQMC method gives accurate results over 
a wide range of physical parameters values. 
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FIG. 1. The spin-spin correlation function as a function of 
the distance R calculated by the ASQMC and the standard 
AFQMC methods at half-filling. The calculations were made 
on the 4x4 lattice. U/t = 4. Open circles and crosses are data 
obtained by the ASQMC and the standard AFQMC methods, 
respectively. 



FIG. 2. The sign expectation value (Sign) as a function of 
At calculated with the ASQMC method. We put Tc — and 
U/t — 8 and p = 0.875. r/t = 10. The calculation is made 
on the 4x4 Hubbard model. 
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FIG. 3. Total energies per site of the two dimensional 
Hubbard model with various system size calculated by the 
ASQMC, ED, the CPQMC, and the AFQMC methods, t = 1 
and [/ = 4. We plot the energies as a function of electron 
density p. Open circles, open squares, and open triangles 
denote results obtained by the ASQMC method with 4x4, 
6x6, and 8x8 lattices, respectively. Crosses, asterisks, and 
closed diamonds denote ED, CPQMC, and AFQMC results, 
respectively. 



FIG. 4. Relative errors of total energies per site calculated 
by the CPQMC and the AFQMC methods. They are com- 
pared with the ASQMC results. The physical parameters and 
fillings are the same as Fig. 3. Crosses and closed diamonds 
denote the relative errors of the CPQMC and the AFQMC 
results, respectively. 
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FIG. 5. Total energies per site as a function of U/t calcu- 
lated by the ASQMC, CPQMC and ED methods. The calcu- 
lations were made on the 4x4 lattice and p = 14/16. Closed 
circles, crosses, and open squares denote ASQMC, ED, and 
CPQMC results, respectively. 



FIG. 6. Errors of total energies per site as a function of U /t. 
The physical parameters and fillings are the same as Fig. 5. 
Closed circles and open squares are errors of the ASQMC 
and the CPQMC results, respectively. 
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FIG. 7. The d-wave superconducting correlation function 
as a function of the distance R calculated by the ASQMC and 
ED methods. The calculations were made on the 4x4 lattice. 
U/t = 4 and p — 10/16. Open circles and crosses are data 
obtained by the ASQMC and ED methods, respectively. 
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